
needs::needs(haven, dplyr, tidyr, ggplot2, readr)

data =  read_dta("../Selling_Out_or_Standing_Firm_Replication Files/Selling_Out_or_Standing_Firm_processlevel.dta")

# CBPS functions

# CBPS
cbps_match = function(formula.match, data, iterations = 2000, method = "over", normalize.weights = T) {
  # set treatment name
  treat.cbps = all.vars(formula.match)[1]
  
  # make treatment a factor
  if(!is.factor(data[[treat.cbps]])) {
    data[[treat.cbps]] = factor(data[[treat.cbps]])
  }
  
  # check for treatment number
  if(n_distinct(data[[treat.cbps]]) > 4) {
    data[[treat.cbps]] = as.numeric(data[[treat.cbps]])
  }
  
  # set ATT
  if(n_distinct(data[[treat.cbps]]) < 3) {
    att = 1
  } else {
    att = 0
  }
  
  # get the formulas
  cbps.form = as.formula(paste(treat.cbps, "~", paste(paste0("`", colnames(model.matrix(formula.match, data))[-1], "`"), collapse = " + ")))
  
  # create the data
  cbps.data = bind_cols(data %>% dplyr::select(!!treat.cbps), as.data.frame(model.matrix(update(formula.match, . ~ . - 1), data)))
  
  # run the model
  cbps = CBPS::CBPS(formula = cbps.form, data = cbps.data, ATT = att, method = method, standardize = T, twostep = T, iterations = iterations)

  # make treatment for two level
  if(!is.factor(cbps$y)) {
    cbps$y = data[[treat.cbps]]
  }
  
  # normalize weights
  if(normalize.weights) {
    cbps$weights = cbps$weights * (nrow(cbps.data) / sum(cbps$weights))
  }
  
  # return
  return(cbps)
}

# Calculate absolute difference in standardized means
adsm_cbps = function(cbps) {
  if(!is.factor(cbps$y)) return(NA)
  
  # data to set means for 
  data.means = bind_cols(data_frame(.treatment = cbps$y, .weights = cbps$weights, .prop.score = 1/cbps$weights), as_data_frame(cbps$x[, -1]))
  
  # gather
  data.means = data.means %>% gather(var, value, -.treatment, -.weights, -.prop.score)
  
  # summarize
  data.means.sum = data.means %>% group_by(.treatment, var) %>% 
    summarize(mean = mean(value), w.mean = weighted.mean(value, .weights), sd = sd(data.means$value[data.means$var == first(var)]), mean.std = mean / sd, w.mean.std = w.mean / sd)
  
  # treatment comparisons
  comps = combn(levels(data.means$.treatment), 2) %>% t %>% as_data_frame
  
  # get the absolute difference in standardized means
  data.abs.dif.std.mean = 
    lapply(1:nrow(comps), function(x) {
      data_frame(t1 = as.character(comps[x, 1]), t2 = as.character(comps[x, 2]), comparison = paste(t1, "vs.", t2),
                 var = data.means.sum$var[data.means.sum$.treatment == t1],
                 Original = abs(data.means.sum$mean.std[data.means.sum$.treatment == t1] - data.means.sum$mean.std[data.means.sum$.treatment == t2]),
                 Balanced = abs(data.means.sum$w.mean.std[data.means.sum$.treatment == t1] - data.means.sum$w.mean.std[data.means.sum$.treatment == t2]))
    }) %>% bind_rows
  
  # spread
  data.abs.dif.std.mean = data.abs.dif.std.mean %>%
    gather(type, diff.mean, -t1, -t2, -comparison, -var)
  
  # return
  return(list(balance = data.means.sum, difference = data.abs.dif.std.mean))
}

#Plot differences 
plot.cbps_balance = function(cbps.bal, var.names = NULL) {
  # the cutoff for ABSMD is generally taken as 0.1: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3144483/#s11title
  cbps.bal$difference$cutoff = if_else(cbps.bal$difference$diff.mean > 0.25, 2, if_else(cbps.bal$difference$diff.mean > 0.1, 1, 0))
  
  # set variable names
  if(!is.null(var.names)) {
    # set var names
    cbps.bal$difference$var.name = var.names[match(cbps.bal$difference$var, names(var.names))] 
    
    # set factor
    cbps.bal$difference$var.name = cbps.bal$difference$var.name %>% factor(levels = var.names[var.names %in% cbps.bal$difference$var.name])
  } else {
    # set var names and factor
    cbps.bal$difference$var.name = cbps.bal$difference$var %>% factor
  }
  
  
  cbps.bal$difference$type_new = factor(cbps.bal$difference$type, levels = c("Original", "Balanced"))
  
  cbps.bal$difference$var.name = as.character(cbps.bal$difference$var.name)
  
  # plot the balance by contrast
  gg.contr = ggplot(cbps.bal$difference, aes(x = diff.mean, y = comparison)) + 
    geom_vline(xintercept = 0.1, color = "darkgrey", linetype = "dashed", size = 1) + 
    geom_point(size = 2, aes(color = factor(cutoff))) + scale_color_manual(values = c("black", "black", "black"), guide = F) +
    #theme_minimal() + 
    facet_grid(cols = NULL, rows = vars(type_new)) + labs(x = "Absolute Difference of Standardized Means", y = "")
  
  # plot by variable
  gg.var = ggplot(cbps.bal$difference, aes(x = diff.mean, y = var.name)) + 
    geom_vline(xintercept = 0.1, color = "darkgrey", linetype = "dashed", size = 1) + 
    geom_point(size = 2, aes(color = factor(cutoff))) + scale_color_manual(values = c("black", "black", "black"), guide = F) + 
    scale_y_discrete(limits = rev(levels(cbps.bal$difference$var.name))) + 
    #theme_minimal() + 
    theme_light(base_size = 12) +
    facet_grid(cols = NULL, rows = vars(type_new)) + labs(x = "Absolute Difference of Standardized Means", y = "")
  
  # return
  return(list("contrast" = gg.contr, "variable" = gg.var))
}

########################
#### Ethnic Support ####
########################

base = ( ~ conf_type_eth + rebstrength_cat + intensity_dummy + osv_state_dum + conventional)
match = update (base, ethnic_all_dum ~ . )

# Main DV
f = update(match, public_proportion2 ~ ethnic_all_dum + . )
d = data %>% ungroup %>% dplyr::select(c("p_id", "dyad_id", all.vars(f))) %>% filter(complete.cases(.))

cbps = cbps_match(match, d, method = "exact")
d$weights = cbps$weights
# calculate absolute difference in standardized means for balance figure
model1 = adsm_cbps(cbps)

eth_match =  glm(f, d, family = quasibinomial('logit'), weights = d$weights)
eth_match_robust = lmtest::coeftest(eth_match, vcov.=sandwich::vcovHC(eth_match, type="HC0"))
eth_match_robust

# Inclusive DV
f = update(match, public_proportion2_inc ~ ethnic_all_dum + . )
d = data %>% ungroup %>% dplyr::select(c("p_id", "dyad_id", all.vars(f))) %>% filter(complete.cases(.))

cbps = cbps_match(match, d, method = "exact")
d$weights = cbps$weights

eth_matchalt =  glm(f, d, family = quasibinomial('logit'), weights = d$weights)
eth_matchalt_robust = lmtest::coeftest(eth_matchalt, vcov.=sandwich::vcovHC(eth_matchalt, type="HC0"))
eth_matchalt_robust

############################
#### External Sanctuary ####
############################

base = ( ~ conf_type_eth + rebstrength_cat + intensity_dummy + osv_state_dum + conventional)
match = update (base, bes_sanctuary_max_duringconf ~ . )

# Main DV
f = update(match, public_proportion2 ~ bes_sanctuary_max_duringconf + . )
d = data %>% ungroup %>% dplyr::select(c("p_id", "dyad_id", all.vars(f))) %>% filter(complete.cases(.))

cbps = cbps_match(match, d, method = "exact")
d$weights = cbps$weights
# calculate absolute difference in standardized means for balance figure
model2 = adsm_cbps(cbps)

ext_match =  glm(f, d, family = quasibinomial('logit'), weights = d$weights)
ext_match_robust = lmtest::coeftest(ext_match, vcov.=sandwich::vcovHC(ext_match, type="HC0"))
ext_match_robust

# Inclusive DV 

f = update(match, public_proportion2_inc ~ bes_sanctuary_max_duringconf + . )
d = data %>% ungroup %>% dplyr::select(c("p_id", "dyad_id", all.vars(f))) %>% filter(complete.cases(.))

cbps = cbps_match(match, d, method = "exact")
d$weights = cbps$weights

ext_matchalt =  glm(f, d, family = quasibinomial('logit'), weights = d$weights)
ext_matchalt_robust = lmtest::coeftest(ext_matchalt, vcov.=sandwich::vcovHC(ext_matchalt, type="HC0"))
ext_matchalt_robust


############################
#### Contraband Funding ####
############################

base = ( ~ conf_type_eth + rebstrength_cat + intensity_dummy + osv_state_dum + conventional)
match = update (base, e_contraband_dum2_max_duringconf ~ . )

# Main DV 
f = update(match, public_proportion2 ~ e_contraband_dum2_max_duringconf + .)
d = data %>% ungroup %>% dplyr::select(c("p_id", "dyad_id", all.vars(f))) %>% filter(complete.cases(.))

cbps = cbps_match(match, d, method = "exact")
d$weights = cbps$weights
# calculate absolute difference in standardized means for balance figure
model3 = adsm_cbps(cbps)

extort_match =  glm(f, d, family = quasibinomial('logit'), weights = d$weights)
extort_match_robust = lmtest::coeftest(extort_match, vcov.=sandwich::vcovHC(extort_match, type="HC0"))
extort_match_robust

# Inclusive DV
f = update(match, public_proportion2_inc ~ e_contraband_dum2_max_duringconf + .)
d = data %>% ungroup %>% dplyr::select(c("p_id", "dyad_id", all.vars(f))) %>% filter(complete.cases(.))

cbps = cbps_match(match, d, method = "exact")
d$weights = cbps$weights

extort_matchalt  =  glm(f, d, family = quasibinomial('logit'), weights = d$weights)
extort_matchalt_robust = lmtest::coeftest(extort_matchalt, vcov.=sandwich::vcovHC(extort_matchalt, type="HC0"))
extort_matchalt_robust


### TABLE 2 and TABLE 2A ### 

labels = c(
  "ethnic_all_dum" =  "Ethnic Support", 
  "b.es_sanctuary_max_duringconf" =  "External Sanctuary",
  "e_contraband_dum2_max_duringconf" = "Contraband Funding",
  "conf_type_eth" =  "Conflict Type",
  "rebstrength_cat" =  "Rebel Strength",
  "intensity_dummy" = "Conflict Intensity",
  "osv_state_dum" = "Government OSV",
  "conventional" = "Conventional")


stargazer::stargazer(eth_match_robust, ext_match_robust, extort_match_robust, covariate.labels = labels,  type = "html", out="../Selling_Out_or_Standing_Firm_Replication Files/Table2-Model3-5.txt")
stargazer::stargazer(eth_matchalt_robust, ext_matchalt_robust, extort_matchalt_robust, covariate.labels = labels,  type = "html", out="../Selling_Out_or_Standing_Firm_Replication Files/Table2A-Model3-5.txt")



##################
#### Figure 1 ####
##################

# read provision level data 

data = read_csv("../Selling_Out_or_Standing_Firm_Replication Files/Selling_Out_or_Standing_Firm_provisionlevel_forFigure1.csv")

# filter regulatory provisions and porvisions that does not include any concession to rebel group
data = data %>% filter (concession != "NA" & concession != -1) 

# labels
data$provision_type_label = case_when(data$provision_type == "amnesty" ~ "Amnesty",
                                      data$provision_type == "demobilize" ~ "Demobilization", 
                                      data$provision_type == "disarm" ~ "Disarmament",
                                      data$provision_type == "integration_civil" ~ "Integration into Civil Admin.", 
                                      data$provision_type == "integration_gov" ~ "Integration into Government",
                                      data$provision_type == "integration_mil" ~ "Integration into Military",
                                      data$provision_type == "integration_police" ~ "Integration into Police",
                                      data$provision_type == "paramilitary" ~ "Dissolution of Paramilitary",
                                      data$provision_type == "prisoner_release" ~ "Prisoner Release",
                                      data$provision_type == "reform_civil" ~ "Civil Admin. Reform",
                                      data$provision_type == "reform_econ" ~ "Economic Reform",
                                      data$provision_type == "reform_electoral" ~ "Electoral Reform",
                                      data$provision_type == "reform_judicial" ~ "Judicial Reform",
                                      data$provision_type == "reform_mil" ~ "Military Reform",
                                      data$provision_type == "reform_party" ~ "Political Party Reform",
                                      data$provision_type == "reform_police" ~ "Police Reform",
                                      data$provision_type == "reform_political" ~ "Political Reform",
                                      data$provision_type == "reform_terr" ~ "Territorial Reform",
                                      data$provision_type == "reintegration" ~ "Reintegration",
                                      data$provision_type == "security" ~ "Security Guarantees",
                                      data$provision_type == "share_econ" ~ "Economic Power Sharing" ,
                                      data$provision_type == "share_terr" ~ "Territorial Power Sharing",
                                      data$provision_type == "transitional_inst" ~ "Transitional Institutions",
                                      data$provision_type == "withdraw" ~ "Withdrawal of Foreign Forces")

# Main DV 
data1 = data
data1$benefit = case_when(data1$public_pos == 1 & data1$allgroupbenefit == 0 ~ "Only Public Benefit",
                          data1$public_pos == 0 & data1$allgroupbenefit == 1 ~ "Only Group Benefit",
                          data1$public_pos == 1 & data1$allgroupbenefit == 1 ~ "Group and Public Benefit")

data1$benefit = factor(data1$benefit, levels = c("Only Group Benefit",  "Group and Public Benefit", "Only Public Benefit"))

# Inclusive DV 
data2 = data
data2$benefit2 = case_when(data2$public_pos2 == 1 & data2$allgroupbenefitinclusive == 0 ~ "Only Public Benefit",
                          data2$public_pos2 == 0 & data2$allgroupbenefitinclusive == 1 ~ "Only Group Benefit",
                          data2$public_pos2 == 1 & data2$allgroupbenefitinclusive == 1 ~ "Group and Public Benefit")

data2$benefit2 = factor(data2$benefit2, levels = c("Only Group Benefit",  "Group and Public Benefit", "Only Public Benefit"))

# Combine for figure
colnames(data2) = colnames(data1)

data1$name = "Main DV"
data2$name = "Inclusive DV"

datacombined = data1 %>% bind_rows(data2)
datacombined$name = factor(datacombined$name, levels = c("Main DV", "Inclusive DV"))

# Figure 

provisiontypes = datacombined  %>% ggplot(aes(x = provision_type_label, fill = benefit )) + 
  geom_bar(position = "stack", stat = "count") + 
  theme_bw() + labs (x = "", y = "Number of Provisions") + 
  theme(legend.direction = "horizontal", legend.position= "bottom" , legend.background = element_rect(), 
        legend.margin = grid::unit(0, "lines"), legend.key.size = grid::unit(0.5, "cm")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_manual(" ", values = c("Only Group Benefit" = "gray85", "Group and Public Benefit" = "gray50", "Only Public Benefit" = "gray15")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + scale_y_continuous(limits = c(0, 50)) + facet_wrap(~name)

jpeg("../Selling_Out_or_Standing_Firm_Replication Files/Figure1.jpg", width = 6.5, height = 9, units = "in", res = 300)
provisiontypes
dev.off()

###########################
#### Online Supplement ####
###########################

###################
#### Figure 1A ####
###################
#Balance Before and After Covariate Balancing Propensity Score (CBPS) Matching 

# labels ----
labels = c(
  "ethnic_all_dum" =  "Ethnic Support", 
  "b.es_sanctuary_max_duringconf" =  "External Sanctuary",
  "e_contraband_dum2_max_duringconf" = "Contraband Funding",
  "conf_type_eth" =  "Conflict Type",
  "rebstrength_cat" =  "Rebel Strength",
  "intensity_dummy" = "Conflict Intensity",
  "osv_state_dum" = "Government OSV",
  "conventional" = "Conventional")

plot_model1 = plot.cbps_balance(model1, var.names = labels)
plot_model2 = plot.cbps_balance(model2, var.names = labels)
plot_model3 = plot.cbps_balance(model3, var.names = labels)


model1_out = plot_model1$variable + ggtitle("Model 4 - Ethnic Support") +
  theme(plot.title = element_text(face="plain", size=14, hjust=0.5))

model2_out = plot_model2$variable + ggtitle("Model 5 - External Sanctuary") +
  theme(plot.title = element_text(face="plain", size=14, hjust=0.5))

model3_out = plot_model3$variable + ggtitle("Model 6 - Contraband Funding") +
  theme(plot.title = element_text(face="plain", size=14, hjust=0.5))


jpeg("../Selling_Out_or_Standing_Firm_Replication Files/cbps_plots.jpg", width = 14, height = 11, units = "in", res = 300)
gridExtra:: grid.arrange (model1_out, model2_out, model3_out, ncol = 3)
dev.off()

